DIPGEOPR_03

class: inverse,left, middle
background-image: url(data:image/png;base64,#background.png)
background-size: cover

##Módulo 3: Análisis masivo de datos satelitales y Google Earth Engine [GEE]

###Google Earth Engine: Aplicaciones en GEE [Parte II]

José A. Lastra<br> <a href="http://github.com/JoseLastra">
Github: JoseLastra</a><br> 
<a href="mailto:jose.lastra@pucv.cl"> jose.lastra@pucv.cl</a> | <a href="mailto:jose.lastramunoz@wur.nl"> jose.lastramunoz@wur.nl</a><br>

.large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Noviembre, 2024</b>] <br>

---
class: center,middle
background-image: url(data:image/png;base64,#labgrs_logo.png)
background-size: 35%

---
## Contenidos

.pull-left[
- Trabajo con colecciones:
<details>

* Importación
  * Filtros temporales y espaciales
  * Creación de una imagen
  * Exportación de datos

</details>
 
- Funciones en GEE
  * Uso de **map()**

- Reducciones
  * Temporales
  * Espaciales
  
- Extracción series de tiempo

]
  
.pull-right[
<img src="data:image/png;base64,#GEE.gif" width="450px"/>

]

---
## Series de tiempo

- GEE tiene disponible un gran número de datasets con información estandarizada en el tiempo, muy útil para análisis de cambios y dinámica espacio-temporal.

- Por ejemplo:
  * GIMMS NDVI (Jul 1981 – Dec 2013),
  * Hybrid Coordinate Ocean Model (Oct 1992 – actualidad),
  * MOD13Q1/MYD13Q1 Versión 6.1 Vegetation Indices 16-Day L3 Global 250 m (Feb 18, 2000 - actualidad),
  * Landsat 5/7/8/9 Collection 2 (Ene 1984 - actualidad)

- Además de información climática de diferentes modelos
- Toda esta información es manipulable y nos permite hacer diferentes análisis además de extraer información.

---
## Extracción de series de tiempo

- Para este ejercicio, aplicaremos algunas herramientas ya usadas: **(i) uso de map()** y **(ii) agregación espacial y temporal** 
- Cargaremos la colección 2 de Landsat 5, 7, 8 y 9 L2 (Surface reflectance) y calcularemos el **Enhanced Vegetation Index** para todas las imágenes disponibles.

- Para filtrar las colecciones haremos uso de una geometría y consideraremos todas las escenas que tengan 70% o menos de cobertura nubosa.

- **Importante**: se puede considerar un umbral más alto o más bajo dependiendo de la zona de estudio y la cantidad de pérdida de datos.

---
## Preparación de funciones

--
- Lo primero a preparar, son las funciones para trabajar con los datos de la colección 2 y los múltiples sensores:
  * Función para estandarizar nombres de bandas
  * Enmascarado de nubes
  * Re-escalamiento de datos ([**más información**](https://www.usgs.gov/landsat-missions/landsat-collection-2-level-2-science-products))
  * Armonización entre sensores ([**más información**](https://openprairie.sdstate.edu/cgi/viewcontent.cgi?referer=https://scholar.google.com/&httpsredir=1&article=1035&context=gsce_pubs))
  * Cálculo del índice

``` js
/* funciones para seleccionar y renombrar bandas
 Landsat 8 y 9 */
var renameOli = function (img) {
  return img.select(
      ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7', 'QA_PIXEL'],
      ['azul', 'verde', 'rojo', 'nir', 'swir1', 'swir2', 'pixel_qa']);
};
// Landsat 5 y 7
var renameEtm = function (img) {
  return img.select(
      ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7', 'QA_PIXEL'],
      ['azul', 'verde', 'rojo', 'nir', 'swir1', 'swir2', 'pixel_qa']);
};
```

---
## Preparación de funciones

--
- **Importante:** la aplicación de quality assesment sobre nuestros datos, dependerá de la zona, el nivel de rigurosidad que queramos en la limpieza, la cantidad de información que estemos dispuesto a perder, entre muchos otros factores.

--
- Enmascarado de datos

``` js
/* limpieza con pixel QA para todos los sensores */
var maskLandsat = function (img){
  return img.mask(img.select('pixel_qa').remap([
          1, // Fill
          5440, // Clear with lows set L5-7
          21824], // Clear with lows set L8-9
          [1,1,1],// Indicar valores a mantener 1 = not masked
            0)); 
};
```

- Basado en: 
  * [**Science Product Guide L4-7**](https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/media/files/LSDS-1618_Landsat-4-7_C2-L2-ScienceProductGuide-v4.pdf) Versión 4.0 Septiembre 2021
  * [**Science Product Guide L8-9**](https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/media/files/LSDS-1619_Landsat8-9-Collection2-Level2-Science-Product-Guide-v6.pdf) Versión 6.0 Mayo 2024

---
## Preparación de funciones

--
- Re-escalado de datos

``` js
// aplicación de factores de escala
var applyScaleFactors = function (image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  return image.addBands(opticalBands, null, true);
};

```

</script>

--
- Solo trabajaremos con la información óptica, pero el escalado también puede aplicarse a la información termal.

---
## Preparación de funciones

--
- Armonización entre sensores

--
- _"Traducción"_ de reflectancias TM/ETM+ a OLI

``` js
/* Armonización Landsat TM/ETM to OLI 
coeficientes provistos por  Roy et al. (2016) */
var coefficients = {
  itcps: ee.Image.constant([0.0003, 0.0088, 0.0061, 0.0412, 0.0254, 0.0172]),
  slopes: ee.Image.constant([0.8474, 0.8483, 0.9047, 0.8462, 0.8937, 0.9071])
};
var etmToOli = function (img) {
  return img.select(['azul', 'verde', 'rojo', 'nir', 'swir1', 'swir2'])
      .multiply(coefficients.slopes)
      .add(coefficients.itcps)
      .addBands(img.select('pixel_qa'));
};
```

---
## __Notas sobre la armonización__

--
- De acuerdo al USGS los productos de la colección 2 vienen corregidos y armonizados:

>_"Yes, Landsat 8/9 are harmonized to Landsat 4/5/7 products as best they can be. There are still the minor differences in the wavelengths collected, but the data is comparable across the different missions."_ USGS, Feb. 2024 (contacto directo).

- Las diferencias en la respuesta espectral o de los índices puede verse condicionada por aspectos como: (i) Índice espectral seleccionado, (ii) ecosistema en estudio, (ii) Cobertura, entre otros.
  
---
class: center, middle
background-image: url(data:image/png;base64,#harmo_NDVI.png)
background-size: 65%

---
class: center, middle
background-image: url(data:image/png;base64,#harmo_EVI.png)
background-size: 65%

---
## Preparación de funciones

--
- Cálculo de EVI

``` js
//Función para calcular EVI 
var EVI = function(image) {
    var evi = image.expression('2.5 * (nir - rojo) / (nir + 6 * rojo - 7.5 * azul + 1)', {
        'nir': image.select('nir'),
        'rojo': image.select('rojo'),
        'azul': image.select('azul')}).rename('EVI');
        
    return image.addBands(evi);
};
```

---
## Preparación de funciones

--
- Todas las funciones creadas las incoporaremos a: **prepDataOli() y prepDataEtm()**

``` js
// Preparación OLI 
var prepDataOli = function (img) {
  var orig = img;
  img = applyScaleFactors(img);
  img = renameOli(img);
  img = maskLandsat(img);
  img = EVI(img);
  return ee.Image(img.copyProperties(orig, orig.propertyNames()));
};
// Preparación TM/ETM
var prepDataEtm = function (img) {
  var orig = img;
  img = applyScaleFactors(img);
  img = renameEtm(img);
  img = maskLandsat(img);
  img = etmToOli(img);
  img = EVI(img);
  return ee.Image(img.copyProperties(orig, orig.propertyNames()));
};

```

</script>

---
## Aplicando funciones a nuestras colecciones

``` js
// aplicando funciones a colecciones
var L5sr = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2")
.filterBounds(roi)
.filter(ee.Filter.lt('CLOUD_COVER',70))
.map(prepDataEtm);
var L7sr = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2")
.filterBounds(roi)
.filter(ee.Filter.lt('CLOUD_COVER',70))
.map(prepDataEtm);
var L8sr = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
.filterBounds(roi)
.filter(ee.Filter.lt('CLOUD_COVER',70))
.map(prepDataOli);
var L9sr = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2")
.filterBounds(roi)
.filter(ee.Filter.lt('CLOUD_COVER',70))
.map(prepDataOli);

```

</script>

---
## Fusión de colecciones

--
- Con nuestros datos preparados, podemos fusionar las colecciones en una sola **ImageCollection**

--
- Esta fusión de colecciones se puede hacer con la función **merge()**

--
- Considerar la función **sort()** para ordenar temporalmente nuestras imágenes dentro de la nueva colección

``` js
var full = L5sr.merge(L7sr).merge(L8sr).merge(L9sr).sort('system:time_start');
print(full, 'coleccion completa');
```

---
## Visualizando EVI

--
- Añadiremos la información al mapa para poder inspeccionar algunos valores.

``` js
// Visualización
var EVI_params = {
min:0,
max:1, 
palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ], 
  bands: 'EVI'};
  
Map.addLayer(full.first(), EVI_params, 'EVI vis');
```

---
class: middle

---
## Construcción series de tiempo

--
- Crearemos un polígono de nombre **zona1** para la extracción mediana de un área de interés.

--
- Emplearemos la función **ui.Chart.image.series.series()**

``` js
// Gráfico de serie de tiempo.
var eviTimeSeries = ui.Chart.image.series(// gráfico utilizando región de referencia
    full.select('EVI'), zona1, // Coleccion y geometria
    ee.Reducer.median(), 30, 'system:time_start') // agregación de mediana
        .setChartType('ScatterChart')// tipo de gráfico
        .setOptions({// opciones
          title: 'EVI',// título de gráfico
          vAxis: {title: 'EVI median', viewWindow: {min: 0, max: 1}},// eje y
          lineWidth: 1,
          pointSize: 2.5,
          series: {
            0: {color: '#4e00f4'}, 
                    }
                    });
// Display.
print(eviTimeSeries, 'serie de tiempo ejemplo');
```

.footnote[
Más detalles sobre la personalización de gráficos la pueden encontrar [__aquí__](https://developers.google.com/chart/interactive/docs)
]

---
class: middle

---
## Construcción series de tiempo

--
- También podemos sacar series de tiempo para múltiples polígonos usando la función **ui.Chart.image.seriesByRegion()**

--
- Veamos un ejemplo con nuestro asset de nombre **zonas**

``` js
var plotEVI = ui.Chart.image.seriesByRegion(
  full, zonas, ee.Reducer.median(), 'EVI', 30)
  .setChartType('LineChart')
  .setSeriesNames(['Pradera','Nativo en río','Paltos','Agua'])
  .setOptions({
    lineWidth: 1,
    pointSize: 2,
    colors: ['#FA0707', '#087605','#97D000','#148CC4'],
    title: 'EVI median',
    vAxis: {title: 'EVI', viewWindow: {min: 0, max: 1}}
  });
  
print(plotEVI, 'Múltiples polígonos');

```

</script>

---
class: middle

---
## Series no imprimibles

--
- Debemos recordar que si tenemos más de 5000 elementos, no podremos imprimir nuestra serie de tiempo.

--
- Lo mismo puede pasar si contamos con demasiados poligonos que llevan a demasiadas agregaciones concurrentes

--
- Tomemos como ejemplo el producto *Optimum Interpolation Sea Surface Temperature (OISST)* de la NOAA v2.1 [2020]

--
- Este producto dispone de datos de de temperatura superficial del oceáno, anomalías y concentraciones de hielo marino desde 1981 hasta la actualidad.

``` js
/* Carga de colección de datos NOAA SST */
var sst = ee.ImageCollection("NOAA/CDR/OISST/V2_1")
.select('sst')
.map(
function(img){
  var rescale = img.select('sst').multiply(0.01);
  return img.addBands(rescale, null, true);
});

// Visualizacion
var sst_vis = sst.first();
print(sst_vis, 'sst visualizacion');
```

// Visualizacion
var sst_vis = sst.first();
print(sst_vis, 'sst visualizacion');
</script>

---

``` js
// parametros visualizacion
var sstParams = {min: -0.180,
  max: 30.00,
  palette: ['040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003']};
Map.addLayer(sst_vis, sstParams, 'Sea Surface Temperature');
```

---
## Series no imprimibles

--
- Tratemos de usar el mismo modelo empleado para la series de tiempo de EVI

``` js

// Extraer la escala nativa de una banda
var scale = sst_vis.select('sst').projection().nominalScale()

// Gráfico de serie de tiempo.
var sst_series = ui.Chart.image.series(// gráfico utilizando región de referencia
    sst, ts_sst, // Coleccion y geometria
    ee.Reducer.mean(), scale.getInfo(), 'system:time_start') // agregación de mediana
        .setChartType('ScatterChart')// tipo de gráfico
        .setOptions({// opciones
          title: 'SST',// título de gráfico
          vAxis: {title: 'SST [°C]'},// eje y
          lineWidth: 1,
          pointSize: 2.5,
          series: {
            0: {color: '#4e00f4'}, 
                    }
                    });
// Display.
print(sst_series, 'serie de tiempo SST');
```

// Extraer la escala nativa de una banda
var scale = sst_vis.select('sst').projection().nominalScale()

---
class: middle

---
## Series no imprimibles

--
- Podemos crear una función, aplicar reductores usando **reduceRegion()** y exportar a csv.

``` js
/* función simple de extracción de series de tiempo largas */
var createTS = function(img){
  var index_data = img.get('system:index');
  var date = img.get('system:time_start');
  var value = img.reduceRegion(ee.Reducer.mean(), ts_sst, scale.getInfo()).get('sst');
  var ts = ee.Feature(null, {'index_data': index_data, 
                             'date': ee.Date(date).format('YYYY-MM-dd'), 
                             'value': value});
  return ts;
};

// Extraccion
var TS = sst.map(createTS);

// Exportar como CSV
Export.table.toDrive({collection: TS, 
selectors: 'index_data, date, value', 
description:'sst_ejemplo',
folder:'sst_NOAA'});

```

// Extraccion
var TS = sst.map(createTS);

// Exportar como CSV
Export.table.toDrive({collection: TS, 
selectors: 'index_data, date, value', 
description:'sst_ejemplo',
folder:'sst_NOAA'});

</script>

---
## Series no imprimibles

--
- ¿Y si tenemos más de un sitio?

``` js
// Definir múltiples sitios como una capa de geometría
var sites = ee.FeatureCollection([
  ee.Feature(ee.Geometry.Point([-75, -42]), {site_id: 'site_1'}),
  ee.Feature(ee.Geometry.Point([-80, -42]), {site_id: 'site_2'}),
  ee.Feature(ee.Geometry.Point([-85, -42]), {site_id: 'site_3'})
]);
```

--
- Aquí creamos un **FeatureCollection** para poder agregar columnas con información

--
- Esto es simil a subir nuestro propio shapefile y usarlo para iterar sobre las colecciones.

---

``` js
// Función para crear series de tiempo por sitio
var createTS = function(site) {
  var site_id = site.get('site_id'); // Identificador del sitio
  // Función interna de reducción y configuracion de la serie
  var tserie = function(img) {
    var index_data = img.get('system:index');
    var date = img.get('system:time_start');
    
    var value = img.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: site.geometry(),
      scale: scale.getInfo(),
      maxPixels: 1e13
    }).get('sst'); // Reducir SST en el sitio
    return ee.Feature(null, {
      'site_id': site_id,
      'index_data': index_data, 
      'date': ee.Date(date).format('YYYY-MM-dd'),
      'value': value
    });
  };
  return sst.map(tserie); // Generar una colección de series de tiempo por sitio
};
```

---
--

``` js
// Aplicar la función a todos los sitios
var TS_all = sites.map(function(site) {
  return createTS(site);
}).flatten(); // Combinar todas las series en una colección

// Exportar como CSV
Export.table.toDrive({
  collection: TS_all,
  selectors: ['site_id', 'index_data', 'date', 'value'],
  description: 'sst_timeseries_all_sites',
  folder: 'sst_NOAA'
});
```

// Exportar como CSV
Export.table.toDrive({
  collection: TS_all,
  selectors: ['site_id', 'index_data', 'date', 'value'],
  description: 'sst_timeseries_all_sites',
  folder: 'sst_NOAA'
});
</script>

---
class: center, middle
background-image: url(data:image/png;base64,#img26.png)
background-size: 65%

---
## Bibliografía complementaria

- Earth Engine Code Editor | Google Earth Engine |. (s.f.). Google Developers. https://developers.google.com/earth-engine/guides/playground

- Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D., & Moore, R. (2017). Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote sensing of Environment, 202, 18-27. [Ver](https://www.sciencedirect.com/science/article/pii/S0034425717302900)

- Mutanga, O., & Kumar, L. (2019). Google earth engine applications. Remote Sensing, 11(5), 591. [Ver](https://www.mdpi.com/2072-4292/11/5/591/htm)

- Zhao, Q., Yu, L., Li, X., Peng, D., Zhang, Y., & Gong, P. (2021). Progress and trends in the application of Google Earth and Google Earth Engine. Remote Sensing, 13(18), 3778. [Ver](https://www.mdpi.com/2072-4292/13/18/3778)

---
class: middle
![logo](data:image/png;base64,#LABGRS-PUCV.png)